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ADAPTIVE ENERGY MINIMISATION FOR /ip^FINITE 
ELEMENT METHODS 

PAUL HOUSTON AND THOMAS P. WIHLER 


Abstract. This article is concerned with the numerical solution of convex 
variational problems. More precisely, we develop an iterative minimisation 
technique which allows for the successive enrichment of an underlying discrete 
approximation space in an adaptive manner. Specifically, we outline a new 
approach in the context of hp—adaptive finite element methods employed for 
the efficient numerical solution of linear and nonlinear second-order boundary 
value problems. Numerical experiments are presented which highlight the 
practical performance of this new hp—refinement technique for both one— and 
two-dimensional problems. 


1. Introduction 


Over the last few decades, tremendous progress has been made on both the math¬ 
ematical analysis and practical application of hnite element methods to a wide range 
of problems of industrial importance. In particular, significant contributions have 
been made in the area of a posteriori error estimation and automatic mesh adapta¬ 
tion. For recent surveys and historical background, we refer to [^15|18||30i32| , and 
the references cited therein. Here, adaptive methods seek to automatically enrich 
the underlying finite element space, from which the numerical solution is sought, 
in order to compute efficient and reliable numerical approximations. The standard 
approach used within much of the literature is to simply undertake local isotropic 
rehnement of the elements (h-refinement). However, in recent years, so-called hp- 
adaptive hnite element methods have been devised, whereby both local subdivision 
of the elements and local polynomial-degree-variation (p-rehnement) are employed. 
These ideas date back to the work by Babuska and co-workers (cf. [4]-[^); see also 
the recent books 11 21 28 29 , and the references cited therein. The exploitation 


of general hp-rehnement strategies can produce remarkably efficient methods with 
high algebraic or even exponential rates of convergence. Moreover, such approaches 
can also be combined with anisotropic rehnement techniques in order to efficiently 
approximate problems with sharp transition features. These techniques enable the 
user to perform accurate and reliable computational simulations without excessive 
computing resources, and with the conhdence that complex local features of the 
underlying solution are accurately captured. 

Many physical processes can be modelled by locating critical points of a given 
(in our setting, convex) energy functional, over an admissible space of functions; a 
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typical example includes quasilinear partial differential equations (PDEs). In this 
article we consider the application of adaptive finite element techniques, employ¬ 
ing a combination of /ip-mesh refinement, to problems of this type. In particular, 
we consider a new and widely applicable paradigm for adaptive mesh generation 
within which we directly seek to construct the hp-fmite element space in order to 
approximate the critical point of the underlying energy functional associated to the 
problem of interest. The simplest such example is the one-dimensional Poisson 
equation on the interval (0,1), subject to a load /; in this case, we seek to minimise 
E(u) = 1/2 fg dx—fg fu dx over an appropriate solution space V (which naturally 
incorporates the boundary conditions). The corresponding standard Galerkin finite 
element approximation of this problem automatically inherits the same energy min¬ 
imisation property with respect to the underlying finite element space I 4 , C V, i.e., 
the finite element solution is the unique minimiser of E(-) over Vh- With this idea 
in mind, a natural approach is to adaptively modify the finite element space I 4 in a 
manner which seeks to directly decrease the energy E, i.e., denoting the new finite 
element solution and finite element space by u'^ and respectively, we require 
that E(uJj) < E{ufi)- By considering an appropriately defined elementwise energy 
functional E),, with k denoting the current element in the underlying computational 
mesh, we devise a competitive refinement strategy which marks elements for refine¬ 
ment. More precisely, in the context of our one-dimensional example, consider an 
/i-refinement strategy which subdivides each element into two sub-elements. We 
may then compute the numerical solution to a local finite element problem posed 
on a local patch of elements which includes the two sub-elements. On the basis of 
this local reference approximation, we may then determine the predicted (elemen¬ 
tal) energy loss if the proposed refinement (i.e., a bisection of each element into two 
sub-elements) is undertaken. Once the predicted energy loss has been computed 
for all elements in the mesh, a percentage of elements with the largest predicted 
energy loss may be identified and subsequently refined. This idea naturally extends 
to the p-refinement setting, whereby, additional higher-order modes are used to 
locally enrich the finite element solution elementwise. With this in mind, we pro¬ 
pose a competitive /ip-adaptive refinement strategy which computes the maximal 
predicted energy loss on each element based on comparing a p-refinement of each 
element (i.e., an isotropic increase of the elemental polynomial degrees by 1) with 
a collection of /i-refinements of the same element (featuring different local poly¬ 
nomial degree distributions), which are selected so as to lead to the same increase 
in the number of degrees of freedom associated with the current element as the 
p-enrichment, cf. [^[^,27|^. A key aspect of this algorithm is the computa¬ 
tion of a local elementwise/patchwise reference solution needed for the definition 
of E),. In order to illustrate the key ideas, for the purposes of this article, we re¬ 
strict our discussion to convex optimisation problems posed on a computational 
domain il C d > 1. However, we point out that this strategy is completely 
general in the sense that it can be applied to any physical problem which may 
be modelled as a critical point of a given energy functional E (including saddle 
point problems, see, e.g., [^). In particular, one of the key advantages of our 
proposed approach is that it naturally facilitates the use of /ip-mesh adaptation, 
and indeed even anisotropic /ip-mesh refinement. By considering an enrichment 
of the finite element space locally using any combination of isotropic/anisotropic 
/i-/p-refinement, an element k can be refined according to the refinement which 
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leads to the maximal predicted energy loss. This is in contrast to standard adap¬ 
tive techniques, whereby elements are marked for refinement according to the size 
of a local a posteriori error indicator. Indeed, in this latter setting, such indicators 
rarely contain information concerning how the local finite element space should be 
enriched, but only indicate that a refinement should be performed. Thereby, alter¬ 
native numerical techniques must be devised which are capable of determining the 
direction of refinement (for anisotropic refinement) or the type of refinement {h- or 
p -). For the latter case, such strategies include regularity estimation [l4|[T^[22p^ , 
use of a priori knowledge M31| , and the computation of reference solutions within 


competitive refinement strategies [3 

11 

13 

16 

17 

22 
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29 , for example. This 

latter class of methods, cf. in particular 11 : 

.2 

>7 

29 and 16 

17 , are very much 


in the spirit of the proposed competitive refinement algorithm developed in this 
article. Finally, we refer to 24 for an extensive review and comparison of many of 
the /ip-adaptive refinement techniques proposed within the literature. 

This article is structured as follows. In Section we briefly present an abstract 
framework for variational problems, and consider an application to quasilinear par¬ 
tial differential equations. Subsequently, in Section the hp-version finite element 
discretisation of such problems is presented, and a new /ip-adaptivity approach is 
developed in detail. The theory will be illustrated with a number of numerical 
experiments on linear and quasilinear boundary value problems in Section Fi¬ 
nally, in Section we summarise the work presented in this article and draw some 
conclusions. 

Throughout this article, we let LP{D), p S [l,oo], be the standard Lebesgue 
space on some bounded domain Z?, with boundary dD, equipped with the norm 
II • IIlp(d). Furthermore, for fc G N, we write W'^’P{D) to signify the Sobolev space of 
order k, endowed with the norm || • || w''.p(d) and seminorm | • \w'^-p(d)- For P = 2 we 
write H^{D) in lieu of moreover, Hq{D) denotes the subspace of H^{D) 

of functions with zero trace on dD. 


2. Variational Problems 


In this section we outline an abstract framework for variational problems in 
Banach spaces, and consider an application to quasilinear boundary value problems. 


2.1. Abstract Minimisation Problem. On a real reflexive Banach space X let 
us consider the minimisation problem 


(2.1) min E(u) = min {F(u) — (?, m)} . 

u£X uGX 

Here, (•,•) is the duality product on A x A', where A' signifies the dual space 
of A, and Z G A' is given. Furthermore, throughout this manuscript, we suppose 
that F : A —>■ K is a continuous and strictly convex functional on A, i.e., 

F(i;i +t{v2 - vi)) < F(z;i) -Ft(F(w2) - F(z;i)) Vt G [0,1] yvi,V2 G A. 

In addition, we make the assumption that F satisfies the coercivity type condition 


( 2 . 2 ) 

where || ■ 

see, e.g., _ 

written in weak form as 


F(u) — {l,u) —>■ -Foo as ||m||js: —>■ oo, 

X is a norm on A. Then, (2.1) possesses a unique minimiser u* G A; 
34 Corollary 42.14]. Furthermore, the problem of finding m* G A can be 


(F'(u*),u) = (Z,i;) ViiGA, 
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provided that F has a sufficiently regular Gateaux derivative F'. 

2.2. Model Problem. We now consider a specific application of the above ab¬ 
stract setting. To this end, let ft C d > 1, be an open bounded Lipschitz 
domain with boundary dfl. We consider the following quasilinear partial differen¬ 
tial equation: 

-V-(/r'(Vw*))+ 5 'K) = /, inff, 

(2.d) 


u = 0 , 


on dfl. 


Here, / = f{x), fj, = /i(Vu), and g = g{u) are given functions, and u* = u*(x) is 
the unknown analytical solution. We suppose that / G for some q> 1. The 

corresponding variational problem reads: 


(2.4) 


min E(m) := min / {/4(Vu) -F g{u) — fu} dx, 

u^X u^X Jq 


where X = for some suitable p > 1. 

With this notation, the following proposition holds. 


Proposition 2.5. Let /i and g from (2.41 be strictly convex and convex, respectively, 
and both continuous on and M, respectively. Furthermore, suppose that, for some 
constants Ci,C 2 > 0, p > 1, and ci, C 2 G K the lower bounds hold, 

(2.6) 

(2.7) giv)>ciV + C 2 , 
as well as the growth conditions 

(2.8) MO<c^2(i + iin, 

(2.9) giv) <C 2 il + \v\P), 

for any ^ S and any rj gM.. Then, for any given f G L^{fl), where Yp-F Y? = 1> 
(2.4) has a unique solution in X = Wq’^(H) as well as on any linear subspace of X. 

Proof. Let us define 

u 1 —^ F(u) := / {g{Xu) + g(u) — fu} da;. 

Jn 

Then, we can cast (2.4) into the abstract framework of ( |2.1[ ), with I = 0 in X'. We 
check the conditions from Section [2T1 separately. To this end, we follow the proof 


presented in 34 Example 42.15]. 

(i) Continuity of F.- Let us consider a sequence {un} C Wg^’^(H), with a limit u G 


lY]’^(r2), i.e., Un ^ u as n 


C Lie) is continuous from [LP{fl)]‘^ to L^(H); see 
Therefore, 


oo. Then, with ( |2.8[ ), the Nemyckii operator 

Proposition 26.6]. 


35 


{/i(VM) — g{Xun)} da; 


< lY(Vu) - M(VM„)||Li(n) -t 0 


as n —)■ oo. Similarly, using (|2.9[), as n —)■ oo, we have that 


0 . 


{giu) - giun)}(ix 


The continuity of u f^fudx follows from / G and from Holder’s 

inequality. Thus, F is continuous. 
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(ii) Strict convexity of F; This simply follows from the strict convexity of /r and 
the convexity of g. 


(hi) Coercivity: According to (2.6) and (2.7), we find that 


F(u) > Cl - fu} dx 

> + J {ciu +C2-fu}dx 

> C'i|w|^i.P(n) “ |ci|||w||Li(n) + C2|f2| - fudx 

where |n| signifies the volume of f2. Employing the Poincare-Friedrich’s and 
Holder’s inequalities, we arrive at 


F(w) > C'||u||^i,p(j^) - (ci + ||/||L<,(n))||u||Lp(n) - C2 
> C'll'“llvvi,p(n) “ (^1 + ll/IU'>(f2))ll^llw"i'!>(a) -^ 2 , 

for some constants ci, C 2 > 0 depending on H. Therefore, it follows that 
E(u) = F(tt) — {l,u) = F(m) —>■ 00 , 


with ||M||wi.p(n) 00 . This is the coercivity condition (2.2). 


The result now follows from 34 Corollary 42.14]. 


□ 


3. /ip-FiNiTE Element Discretisation 

Consider now a linear subspace A„ C X with dim(X„) = n < 00 . Then, by 
our previous assumptions on F, solving the finite dimensional convex optimisation 
problem min„gx„ E(m) for the unique minimiser m* G A„ results in an approxi¬ 
mation oi u* G X from (2.1) with m* « u*. This is the well-known Ritz method. 


Equivalently, in weak form, we may seek u* G Xn such that the Galerkin formula¬ 
tion 

{F\u*J,v) = {l,v) VuGA„ 


is satisfied; cf. 34, § 42.5]. 


For the purposes of discretising our model problem (2.3), we will focus on an 


hp-Umte element approach. To this end, let us first introduce some notation: We 
let T = {«:} be a subdivision of the computational domain H into disjoint open 
simplices such that H = lj„g 7 -k and denote by /i„ the diameter of k G T; i.e., = 

diam(K). In addition, to each element k G T we associate a polynomial degree 
Pk. > I, and collect the in the polynomial degree vector p = [p,^ : k G T]- With 
this notation we define the /ip-finite element spaces by 


'^{T,p) = {v G : v\kG 

Vo(r,p) = V(r,p)niJi(E!), 


Xk) vkg r}. 


where, for p > 1, we denote by Pp(k) the space of polynomials of total degree p 
on K. 


The /ip-version finite element approximation of the variational formulation (2.1) 
is given by: Find the numerical approximation G Vo(T, p) such that 

E«p) = min E(u), 

^ «GVo(r,p) 
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where E is defined in (2.4), or equivalently, provided that the (weak) derivatives 
and g' belong to in weak form: Find S '^o{'T,p) such that 


(3.1) 

Here, 


an (uhp ,v) = iniv) Vu e Vo(T, p)■ 


an{w,v) := 


£q{v) := 


{p'{Vw) ■ Vv + g'{w)v} dx, 
fvdx, 


Jn 


This is the hp-finite element discretisation of (2.3). 


w,v £ V{T,p), 
vGV{T,p). 


4. lip-ADAPTIVITY 


The goal of this section is to design a procedure that generates sequences of hp- 
adaptively refined finite element spaces in such a manner as to minimise the error 
in the computed energy functional E. In the context of hp-version finite element 
methods, the local finite element space on a given element k, k £ T, may be enriched 
in a number of ways. In particular, traditional hp-adaptive finite element methods 
typically make a choice between either: 


• p-refinement: The local polynomial degree p^ on k is increased by a given 
increment, pinc^ Pk. ^ Pk + Vine- Typically, a value of pi^c = 1 is selected. 

• /i-refinement: The element k is divided into a set of new sub-elements, 
such that K = ljr=i Here, will depend on both the type of element to 
be refined, and the type of refinement employed, i.e., isotropic/anisotropic. 
For isotropic refinement of a triangular element k in two-dimensions, we 
have Uk = 4. The polynomial degree may then be inherited from the parent 
element k, i.e., we set p^. = p^, for i = 1,..., n^. 

Motivated by the work presented in 


27 , cf. also 1111112, 29 , for example, in this 


article, we consider a competitive refinement strategy, whereby on each element k 
in T, we estimate the predicted reduction in the local contribution to the energy 
functional E based on either employing p-refinement, with pi^c = 1, together with 
a series of hp-refinements, which lead to the same number of degrees of freedom 
as the p-enrichment. In contrast to standard /i-refinement, where the subdivided 
elements inherit the polynomial degree of their parent, cf. above, in this latter case, 
the distribution of the polynomial degrees on the resulting sub-elements is possibly 
non-uniform. 


4.1. Motivation. The key to the forthcoming hp-refinement strategy is to esti¬ 
mate the predicted reduction in the energy functional locally on each element in 
the finite element mesh T. With this in mind, we must first rewrite E as the sum 
of local contributions on 'T. Given that E is simply defined as an integral over O, 
then clearly, we may write 

where is defined in an analogous fashion to E, with the integrals over Q being 
restricted to integrals over k, k G T- 

However, while the above definition of the local energy functionals E„ seems 
entirely natural, there is no guarantee that the computed error will converge op¬ 
timally based on locally minimising E^ over each k in T. In order to investigate 
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this issue further and to motivate the idea proposed in this article, let us consider 
the following second-order linear self-adjoint partial differential equation: Find u* 
such that 


—Au* + u* = f in n, u* = 0 on dfl. 


Thereby, we have that fi{\7u) = 1 / 2 !Vup and g{u) = ^ and X = (i.e., 

p = q = 2). In this setting, the (global) energy functional from (2.41 may be written 
in the form 

E(m) = T^aniu, u) - £n{u). 


Moreover, we may define the associated energy norm: Ijuljg := an{u,u). Given 
the energy norm, exploiting the symmetry of the bilinear form aa(-, •), we immedi¬ 
ately deduce the following relationship between the error in the computed energy 
functional E, and the error measured in the terms of the energy norm Ij-Hg, namely: 


E(u*) - E«p) 



2 

E ■ 


Thereby, on a global level, reduction of the error in the energy functional E naturally 
leads to a reduction in the energy norm of the error. 

In order to repeat this argument on a subset V d fl, we now suppose that the 
boundary datum g is given and seek u* G H^(T)) such that u*\gx> = 9 and 


(4.1) 


aviu* ,v) = £v{v) \IvGHI[V). 


Here, ax){-,-) and iv{') are defined in an analogous manner to an(-,-) and 
respectively, with the domain of integration restricted to V. In this case, writing Tij 
and px> to denote the finite element sub-mesh and polynomial degree distribution 
over V, respectively, the finite element approximation is given by: Find ujjp € 
V{Tv,Pv) such that Uhpla'D = and 

a-Diu^-pyv) = £t>{v) yv G Vo{Tt>,Pv), 

where lig denotes a piecewise polynomial approximation in iJ {dV) of the Dirich- 
let datum g. Thereby, writing E-p to denote the restriction of the energy functional 
E over V, i.e., 

Ep(u) := ^aTi{u,u) - £t>{u), 

we deduce the following identity: 

Ep(u*) - Ep(Mhp) = -^aviu* - ulp,u* - u^p) + av{u*, u* - u^p) - ^p(u* - Uhp)- 

Employing integration by parts we deduce that 

r Qn* 

ap(u*,u*-<p)-^p(M*-<p) = / p-(u*-<p)ds 

Jdv c'np 

= / ^'(Vu*) • np(u* - ■Uhp)ds, 

Jdv 

where np denotes the unit outward normal vector on the boundary dT> of the 
domain V. Thereby, 

= -l^aT,{u*-u^p,u*-ulp)+ [ Ai'(VM*)-np(u*-<p)ds. 

^ JdV 


(4.2) Ep(u*)-Ep«p) 
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Stimulated by (4.2), we define the local energy functional E'p(-) by 

(4.3) Eti{v) ■■= E.xi{v) - [ ■ nxivds. 

Jdv 

With this definition, we immediately deduce the following relationship between the 

error in the local energy functional and the error measured in terms of the local 

energy norm, namely, 

~ 111 112 

Ex.(u*) - Ex,«p) = -- ||u* - Uhp||E_.p , 

where, for w € we let HwUlx) aTi{w,w). Moreover, if we consider the 

evaluation of the above local energy functional on each element n, hi G 'T, then we 
note the following consistency condition holds 

E(u) = ^ E.(n). 


Let us now write 


S V{Tt>,Pv) and G V{T-[,,p'x>) to denote two finite el¬ 
ement approximations to ( |4.1[ ) based on employing the computational meshes T® 
and 7^, respectively, with polynomial degree vectors px> and p'j,, respectively. As¬ 
suming the finite element space V{T^,p'-jy) represents an enrichment of the original 
one V(7 i5,Px)), we deduce that the expected reduction in the error in the energy 
functional dehned over "D satisfies the equality 

E^«p) - - (E^(U*) - Ep«p)) 

1 II * l|2 1 * */ 2 

= 2^ -«hp||E,^-2 - -%p 

Hence, by employing the modihed local definition of the energy functional E-o de¬ 
fined over the subdomain I?, we observe that the expected reduction in Ex> is directly 


related to the reduction in the energy norm of the error over V. The equality (4.4) 
will form the basis of the proceeding /ip-adaptive refinement algorithm. 

4.2. Competitive hp—refinement strategy. In this section, we develop an hp- 
adaptivity algorithm based on employing a competitive refinement strategy on each 
element n in the computational mesh T. The essential idea is to compute the max¬ 
imal predicted energy reduction EK(uj)p) — E„(m* on each element k G T, where 
M^p is the (global) finite element element solution defined by (3.1), and is the 
(local) finite element approximation to the analytical solution u* evaluated on a 
local patch of elements neighbouring k, subject to a given p-//ip-refinement. Em¬ 
ploying the forthcoming notation u* loc will either represent u* p S V{Ti^,Pp), cf. 


(1^, or It* e V(T;)^ef,Php,), i = 1,.. ■,N^,hp, cf. (|T^, corresponding to either 


a local p- or /ip-refinement of element k, respectively. Elements with the largest 
maximal predicted decrease in the local energy functional are then appropriately 
refined. However, before we proceed, we first note that the boundary correction 
term included within the definition of the local energy functional Ek(-)i cf. (4.3) 


with D replaced by k, is not computable since it directly assumes knowledge of the 
unknown analytical solution u*. With this in mind, we replace u* by an approxi¬ 
mate reference solution, cf. [Tlj[^ . However, in contrast to these citations, for the 
purposes of the current article we simply compute local reference solutions, rather 
than global ones. More precisely, given k G T, we first construct the local mesh 
comprising of k and its immediate face-wise neighbours, cf. Fig. [^b). Given 
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(a) 



Figure 1. Local element patches in two-dimensions, when trian¬ 
gular elements are employed, (a) Original element n (assumed to 
be an interior element); (b) Mesh patch tA, which consists of the 
element n and its neighbours; (c) Mesh patch which is con¬ 

structed based on isotropically refining k (red refinement) and on 
a green refinement of its neighbours. 



we then uniformly (red) refine element n into Uk sub-elements; the introduc¬ 
tion of any hanging nodes may then be removed by introducing additional (green) 
refinements, or alternatively, by simply uniformly refining all elements in the sub¬ 
mesh ■ For the purposes of the article, in two-dimensions, we exploit the former 
strategy, purely on the basis of reducing the number of degrees of freedom in the 
underlying local finite element space. Denoting the resulting finite element mesh 
by cf. Fig.^c), we construct the finite element space V(7)(^ef, Pref), where 

Pref U' = Pk + 1 for all k' e 7;^ef Writing V{k) = the elementwise 

reference solution may be computed as follows: Find u* j-ef ^ V(7Aef:Pref) such 
that la-D(K) = <plax>(«;) and 

(4.5) = ^'D(k){v) Vu e Vo(7;^ef.Pref)- 

On the basis of the computed reference solution, we define the approximate local 
energy functional on n, k G T, as follows: 

(4.6) E(,(u) := E^(u) - / )} • ds, 

J Bk, 

where denotes the unit outward normal vector to the boundary da of k, and 
denotes the average operator. More precisely, given two neighbouring elements k'^ 
and K ~, let x be an arbitrary point on the interior face given by F = dn'^ H dK~. 
Given a vector-valued function q which is smooth inside each element , we write 
to denote the traces of q on F taken from within the interior of , respectively. 
Then, the average of q at a: S F is given by -g^qj = ^(q'*" + q~). On a boundary 
face F c 90, we set -gfqj- = q*^. 
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Figure 2. Polynomial degree distribution employed for the com¬ 
petitive hp-refinements: (a) One-dimension; (b) Two-dimensional 
triangular element. 


With the definition of E^(-) given in (4.6), we now outline the proposed compet¬ 
itive refinement strategy on element k, k G T- Firstly, we compute the predicted 
energy functional reduction when p-refinement is employed, i.e., 




(4.7) AE',,p := EU<p) - E',«,p), 

where u* p is the solution of the local finite element problem: Find u* p 
such that ■u*,p|aD(„) = ■Uhpbx)(„) and 

(4.8) Vu e Vo(7;^,Pp); 
here, = p^ -I- 1 for all k' G . 

Secondly, we also consider a sequence of competitive hp-refinements, such that 
the number of degrees of freedom associated with the finite element space defined 
over K is identical to the case when pure p-refinement has been employed. Here, 
for each element k S T, we again exploit the same local mesh employed 

for the computation of the local reference solution u* . Then for the elements 
which result from the isotropic refinement of k, we employ local polynomial degrees 
p„., i = for the remaining elements stemming from the refinement of 

the neighbours of k, we simply set the local polynomial degree equal to p^, cf. 
Fig. i For example, in one-dimension, following [11[|29| , given an element k with 
polynomial degree p^, an enrichment of p^ —?► p^ + 1 gives rise to p^ + 2 degrees 
of freedom associated with k. On the other hand, we can now consider the case 
when K is uniformly subdivided into two sub-elements ki and K 2 , i.e., = 2, with 

associated polynomial degrees p^^ and p^^, respectively. To ensure that the number 
of degrees of freedom in the underlying /ip-refined finite element space defined over 
Ki and K 2 is identical to the case when pure p-enrichment is undertaken, we require 
that 


Pki +Pk2 =Pk + 1. 

Hence, there are = Pk, /ip-competitive refinements and one p-refinement in 

one-dimension. 

In higher-dimensions, the construction of the competitive /ip-refinements is un¬ 
dertaken in an analogous manner. For simplicity, we focus on the two-dimensional 
case when triangular elements are employed. Then for the elements which re¬ 
sult from the isotropic refinement of n, we employ local polynomial degrees p^^, 
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Figure 3. Number of competitive /ip-refinements, versus 

the local polynomial degree when a triangular element k is 
isotropically refined. 


i = 1,..., = 4; as before, the local polynomial degree of the remaining elements 

stemming from the refinement of the neighbours of k is set equal to Let us 
signify the set of all such polynomial degree distributions on T^ef by ^k,p^ ■ Given 
that the full space of polynomials has been employed for the p-refinement, the 
number of degrees of freedom associated with k is + 2)(pk + 3). Then, for 

an arbitrary polynomial degree distribution {pKi}i-i for the sub-elements 
of K, the number of degrees of freedom associated with k is 

3 ^ 4 

6 + ^ [min(p«,,p«J -1 + 2(p«, - 1)] + a - 2), 

where we have assumed that K 4 is the sub-element located at the interior of cf. 
Fig. mb). Thereby, we select the set of /ip-refinements which satisfy the condition 

3 ^4 ^ 

6-f^ [min(p^,,p^J - 1-b 2(p^, - 1)]-b - ^(p«, - l)(p,^. -2) = -(pK-b2)(p^-b3). 

z=l ^ i=l 


Analogous expressions can also be determined for different element types, other 
kinds of refinement, e.g., anisotropic refinement, as well as in higher-dimensions. 
The precise number of competitive /ip-refinements, denoted by is not possible 

to determine in a simple closed form expression; instead, iV^^hp can be precomputed 
for any polynomial order. To this end, in Fig. we present the number of combina¬ 
tions of local polynomial degrees {pKi}i=i with respect to p^ in the above setting, 
i.e., for the case of isotropic rehnement of a triangular element in two-dimensions. 
We notice that the number A^K,hp of possible p-configurations is, not surprisingly, 
growing as p^ increases. In view of this observation we remark that, although 
the subsequent local discrete problems defined on each corresponding (patchwise) 
hp-finite element space, cf. (4.9) below, are extremely inexpensive to compute, 
and moreover are trivially parallelisable, from a practical point of view, it might 
be computationally beneficial to limit the number of samples to a certain preset 
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maximum iVmax- For example, a random selection of A^max samples may be con¬ 
sidered, cf. Section below; alternatively, a more sophisticated strategy selecting 
polynomial degree distributions with limited variations could be employed. 

We now write V(7^eoPhp.), i = 1,..., iVK,hp, to denote the finite element space 
based on employing the local (refined) mesh and some local polynomial degree 

distribution php^ G ■ Thereby, the following competitive /ip-refinements may 
be defined: Find G V,Php.) such that <_hpjax>(«;) = KplorJin) and 

(4.9) a'D(K)«,hp.>'(^) = ^ v ( k ){ v ) Vp G Vo(7;^ef>Php,), 

for i = 1,..., Nf^^hp- For each local competitive /ip-refinement, we compute the 
estimated local energy reduction 

(4T0) hp. := EU<p) - 

for i = 1,... ,Nfillip. In this way, for each element k G T, we may compute the 
maximum local predicted error reduction 

(4.11) = max < AE(, , max AE(, jj > , 

’ i=l,...,Af,i,hp ’ ') 

. Finally, we refine the set of elements k G T which satisfy 

the condition 


with AE' from (4.7 


(4.12) 


AE' 


> 0maxAE' 


k^'T' 


where 0 < 0 < 1 is a given parameter, cf. 11,29. On the basis of |11[|29|, through¬ 


out this article, we set 9 = 1 / 3 . The above competitive /ip-refinement strategy is 
summarised in Algorithmic 


5. Numerical Examples 


In this section we present a series of numerical experiments to demonstrate the 
practical performance of the proposed /ip-adaptive refinement strategy outlined in 
Algorithm [C 


5.1. Example 1: Linear Elliptic Problem. In this first example, we consider 
a one-dimensional problem defined over the domain O = (0,1). Moreover, we set 
^J-{ux) = 1 / 2 e > 0, g{u) = 1 / 211 ^, and f{x) = 1; this is equivalent to solving 
the linear elliptic boundary value problem: 

—su%^ -I- u* = 1, X G 


subject to homogeneous Dirichlet boundary conditions. We note that the analytical 
solution is given by 


^(x) = 


- 1 


G/vT _ o-VvT 


e"/'^ - 1 - 


1 - eF^ 
.VvT _ o-Vu?" 


p ^/e 


-El. 


In particular, for 0 < e <C 1 
the vicinity of x = 0 and x = 1, cf. 33 


the analytical solution u* contains boundary layers in 
we set e = 10“^. 


as in 


33 


In Fig. I^we illustrate the performance of the proposed /ip-adaptive algorithm, 
cf. Algorithmic based on a starting mesh consisting of 4 elements, with the initial 
polynomial degree p = [1,1,1,1]. Here, we have plotted the error in the underlying 
energy functional E, together with the energy norm || • ||e and L^(H) norm of the 
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Algorithm 1 Competitive /ip-adaptive refinement procedure 

1 : Choose a coarse initial mesh To of fl and a corresponding low-order starting 
polynomial degree vector pq. Set n = 0. 


Solve ( phlj ) for eV{Tn,Pn)- 
for each element k G Tn do 

Construct the local reference mesh T^gf • 

Compute the local finite element reference solution m* j-ef G V(7^ef!Pref) 


satisfying (4.51. 


6 : Compute the local finite element p-enriched solution u* p G V(7^,Pp) sat 


isfying (4.8), together with the corresponding predicted energy functional 


reduction AE'^ p, cf. (|4.7|). 


for i = 1,..., iV„,hp do 

Compute the local competitive hp-refined finite element solutions 
^(T^renPhpJ satisfying (4.9), together with their respective 


K.hp- 


predicted energy functional reduction AE'^j^p defined in (4.10). 
end for 

Compute the maximum local predicted error reduction AE^ cf. (4.11). 

end for 

Determine the set of elements /C„ which are flagged for refinement, based on 


the criterion (4.12). 


13: Perform p- or hp-refinement on each k G /C„ according to which refinement 


takes the maximum in (4.11). This results in a refined global finite element 
space V(7;i+i,p„+i). 

14: Set n ^ n -I- 1, and goto Line 2 . 

15: After sufficiently many iterations have been performed output the final solution 
<p G V{Tu,Pu)- 


error, with respect to the total number of degrees of freedom employed within the 
finite element space V(T,p), on a linear-log scale; here, 


= / 

JO 


dx. 


From Fig. I^a), (b), & (c), we observe, that after an initial transient, the conver¬ 
gence lines for each error measure become (on average) straight, thereby indicat¬ 
ing exponential convergence of the quantities |E(w*) — E(u^p)|, ||m* — u((p||E, and 
||u* —M(Lp||i 2 (Q), respectively, as V(T,p) is adaptively enriched. Finally, in Fig. gd) 
we show the hp-mesh distribution after 9 adaptive refinements. Here, we observe 
that the algorithm clearly identihes the location of the boundary layers present in 
the analytical solution u*\ indeed, in these regions, local subdivision of the mesh 
has first been employed, followed by subsequent p-enrichment, cf. 


33 


5.2. Example 2: Strongly monotone quasilinear PDE. In this second exam¬ 
ple, we let H be the L-shaped domain (—1,1)^ \ [0,1) x (—1,0], and set 

p(Vn) = ^(|Vw|2-e-l^“l'). 
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Figure 4. Example 1. Comparison of the error with respect to the 
number of degrees of freedom: (a) |E(u*) — E(u^p)|; (b) ||n*— uJ^p||e; 
(c) ||m* — uj^p||i 2 (f 2 ). (d) hp-Mesh distribution after 9 adaptive 
refinements. 


Thereby, the corresponding Euler-Lagrange equation for the underlying minimisa¬ 
tion problem corresponds to the strongly monotone quasilinear PDE given by: 

(5.1) -V- Vu*) =/, inn. 


We select / and appropriate inhomogeneous Dirichlet boundary conditions so that 


the analytical solution to (5.11 is given by 

V3, 


u = r' sm 


b) 


where (r, p) denote the system of polar coordinates, cf. [lOl[^ , for example. 

Selecting the energy norm || • ||e to be the standard H^{n) norm, in Fig. we 
again present the convergence history of the error in the computed energy func¬ 
tional E, together with ||u* — WepHe, and ||u* — as the finite element 

space is /ip-adaptively rehned. On a linear-log scale (where the horizontal axis 
measures the third root of the total number of the degrees of freedom, cf. [^), we 
again observe exponential rates of convergence, in the sense that asymptotically the 
convergence lines become roughly straight. In addition, in Fig. [^we also present 
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(a) (b) 



Figure 5. Example 2. Comparison of the error with respect to the 
third root of the number of degrees of freedom: (a) | E (m* ) — E (w^p) |; 
(b) ||m*-<pI|e; (c) ||u*-<p||L2(n). 


analogous results in the case when a Monte Carlo (MC) approach is employed 
to limit the number -/Vmax of hp-refinement samples considered on each element. 
More precisely, we randomly select samples based on employing iV^ax = 10 and 
Nmax = 15; in each case two typical realisations are presented. Here, we observe a 
slight degradation of the rate of convergence in each of the above error quantities 
as our hp-refinement procedure progresses, as we would expect; however, in each 
case exponential convergence is retained when this simple selection principle is ex¬ 
ploited. As noted in Section |4.2| more sophisticated selection principles may also 
be employed. 

The final hp-mesh distribution is depicted in Fig. here, we see that the com¬ 
putational mesh has been largely refined in the vicinity of the re-entrant corner 
located at the origin. In addition, we see that the polynomial degrees have been 
increased away from the origin, since the underlying analytical solution is smooth in 
this region. In particular, we observe that the refinement algorithm has generated 
an hp-mesh distribution which is symmetric with respect to the line X 2 = —xi. 
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Figure 6. Example 2. (a) ft.p-Mesh distribution after 18 adaptive 
refinements; (b) Zoom of (a). 


5.3. Example 3: p— Laplacian. In this final example, for p > 1, we consider the 
p-Laplacian problem 

(5.2) - V • (|Vw*|P-2Vu*) =/, in^l = (0,l)^ 


subject to inhomogeneous Dirichlet boundary conditions. We point out that in this 
setting, (5.21 corresponds to the Euler-Lagrange equation for the energy minimi¬ 


sation problem 


min 

MGWo'''(n) 


- j |Vu|Pda:— f fudx 
P Jn Jn 


i.e., we have /i(VM) = g = 0. We select /, and impose suitable 

inhomogeneous Dirichlet boundary conditions, so that the analytical solution of 


(5.2) is given by 


^(x) = ■ 


a > 0. 
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Figure 7. Example 3. Comparison of the error with respect to the 
third root of the number of degrees of freedom: (a) |E(m*) —E(M^p)|; 
(b) |w*-<p|vyi.3(n); (c) ||m* - <pI|l3(o). 


As in throughout this section, we set p = 3 and a = which implies that 
u* G where /3 = is/e and e > 0 is arbitrarily small. 

In Fig.j^we plot |E(u*) - E(M*p)|, ||u* - M*p||i4,i,3(n), and ||u* - <p||L3(n), with 
respect to the third root of the number of degrees of freedom in V{T,p). As in 
the previous examples, we again observe exponential convergence of each of the 
above error measures, as the finite element space is hp-adaptively modified. Here, 
we also consider the case when A^max = 30 random samples are selected; as in 
the previous example, we again see that exponential convergence of each of the 
above error quantities is retained, though the rate of convergence is inferior when 
compared to the case when all potential trial hp-refinements are considered. The 
final hp-mesh distribution is shown in Fig. as in the previous examples, the 
adaptive algorithm clearly identifies the location of the singularity present within 
the analytical solution u*, whereby h-refinement is undertaken. 

6. Conclusions 

In this article, we have proposed a novel hp-adaptive refinement procedure for 
application to the finite element approximation of convex variational problems. 
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0 xk :- 

0 0.2 0.4 0.6 0.8 1 


(a) 



Figure 8 . Example 3. (a) /ip-Mesh distribution after 23 adaptive 
refinements; (b) Zoom of (a). 


In particular, the underlying adaptive algorithm exploits a competitive refinement 
technique which seeks to maximise the decrease in the elemental contribution to the 
total energy based on employing local p- and hp-enrichments of the finite element 
space. Whilst our approach has been successfully applied to a range of second- 
order quasilinear problems in both one- and two-dimensions, we emphasise that it 
is immediately extensible to more general variational-based PDE problems. Future 
work will be concerned with exploiting anisotropic hp-mesh adaptation. 
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